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Abstract. We undertake a systematic, direct numerical simulation (DNS) of the 
two-dimensional, Fourier-truncated, Gross-Pitaevskii equation to study the turbulent 
evolutions of its solutions for a variety of initial conditions and a wide range of 
parameters. We find that the time evolution of this system can be classified into four 
regimes with qualitatively different statistical properties. First, there are transients 
that depend on the initial conditions. In the second regime, power-law scaling regions, 
in the energy and the occupation-number spectra, appear and start to develop; the 
exponents of these power-laws and the extents of the scaling regions change with 
time and depended on the initial condition. In the third regime, the spectra drop 
rapidly for modes with wave numbers k > kc and partial thcrmalization takes place 
for modes with k < kc\ the self-truncation wave- number kdt) depends on the initial 
conditions and it grows either as a power of t or as logi. Finally, in the fourth regime, 
complete-thermalization is achieved and, if we account for finite-size effects carefully, 
correlation functions and spectra are consistent with their nontrivial Berezinskii- 
Kosterlitz-Thouless forms. 



PACS numbers: 47.27.Ak, 47.27.Gs, 47.37.-|-q, 67.25.dk, 67.25.dj 



1. Introduction 



The elucidation of the nature of superfiuid turbulence, which began with the pioneering 
studies of Feynman [1] and of Vinen and Hall [2H6], has continued to engage the attention 
of experimentalists, theoreticians, and numerical simulators [714T3]. Experimental 
systems, in which such turbulence is studied, include the bosonic superfiuid ^He, 
its fermionic counterpart ^He, and Bose-Einstein condensates (BECs) of cold atoms 
in traps and their optical analogues; for representative studies we refer the reader 
to [T1H23]. Theoretical and numerical studies have used a variety of models to 
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study superfluid turbulence; these include the two- fluid model [SUES], Biot-Savart- 
type models with [261 127] or without [281 [29] the local-induction approximation, and 
the Gross-Pit aevskii (GP) or nonlinear Schrodinger (NLS) equations [301 [3T|. These 
models have been studied by a combination of theoretical methods, such as wave- 
turbulence theory [30ti33] . and numerical simulations [3H - H0] . Most of these studies have 
been carried out in three dimensions (3D); numerical simulations of two-dimensional 
(2D) models for superfluid turbulence have been increasing steadily over the past few 
years [HHS]. Here we undertake a systematic direct numerical simulation (DNS) of the 
dissipationless, unforced, Fourier-truncated, 2D, GP equation with a view to identifying 
what, if any, features of the turbulent evolution of the solutions of this equation are 
universal, i.e., they do not depend on initial conditions. Some, though not all, parts of 
our results are contained in earlier simulations [4TH47] . The perspective of our study 
is different from earlier studies of the 2D GP equation; in particular, we elucidate 
in detail the dynamical evolution of this system and examine the various stages of its 
thermalization; in this sense our work is akin to recent studies of thermalization in Euler 
and other hydrodynamical equations [181450] : a similar study for the 3D GP equation 
has been carried out by Krstulovic and Brachet [381 IST] . 

It is useful to begin with a qualitative overview of our principal results. We find 
that the dynamical evolution of the dissipationless, unforced, 2D, Fourier-truncated 
GP equation can be classified, roughly, into the following four regimes, which have 
qualitatively different statistical properties: (1) The first is the region of initial 
transients; this depends on the initial conditions. (2) This is followed by the 
second regime, in which we see the onset of thermalization; here the energy and 
occupation-number spectra begin to show power-law-scaling behaviours, but the power- 
law exponent and the extents of the scaling regions change with time and depend on 
the initial conditions. (3) In the third regime, which we call the region of partial 
thermalization, these spectra show clear, power-law, scaling behaviours, with a power 
that is independent of the initial conditions, and, at large wave vectors, an initial- 
condition-dependent, self-truncation regime, where spectra drop rapidly; (4) finally, in 
the fourth regime, the system thermalizes completely and exhibits correlation functions 
that are consistent with the predictions of the Berezinskii-Kosterlitz-Thouless (BKT) 
theory [ITIESHM], if the simulation domain and simulation time are large enough. 
Although some of these regimes have been seen in some earlier numerical studies of the 
2D GP equation, we are not aware of any study that has systematized the study of these 
four dynamical regimes. In particular, regime 3, which shows partial thermalization 
and self-truncation in spectra, has not been identified in the 2D, Fourier-truncated, GP 
equation, even though its analogue has been investigated in the 3D case [32 l [38 | [5T] . 

The remaining part of this paper is organised as follows. In section [21 we describe 
the 2D, GP equation and the different statistical measures we use to characterize 
turbulence in the Fourier-truncated, 2D, GP equation (section l2.ip : the details of our 
numerical methods and initial conditions are given in section 12.21 In section |3l we 
present our results; these are described in the four subsections 13 . 1113^ that are devoted, 
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respectively, to the following: (a) the temporal evolution of the energy components, 
velocity-component probability distribution functions (PDFs), and the population Nq 
in the zero-wave-number mode; (b) the statistical characterization of the first two 
regimes of the dynamical evolution (by using various energy and the occupation-number 
spectra for different initial conditions); (c) a similar statistical characterization, as in 
subsection 13.21 but for the regime with partial thermalization, and the study of the 
nature of the growth of the self-truncation region; (d) the final, completely thermalized 
state of the Fourier-truncated, 2D, GP equation. Section H] contains our conclusions. A 
note on the units used for the GP equation and the details of some analytical calculations 
are presented in Appendix A and Appendix B, respectively. 



2. Model, Initial Conditions, and Numerical Methods 

In this Section, we describe the 2D, GP equation. We define all the statistical measures 
that we use to characterize the time evolution of this equation, given the three types of 
initial conditions that we describe below. We also describe the numerical methods, and 
computational procedures that we use to solve this equation. 



2.1. The Gross- Pitaevskii Equation 



The GP equation, which describes the dynamical evolution of the wave function ip of a. 
weakly interacting 2D Bose gas at low temperatures, is 



^^^ = -VV(x,t)+^?HV(x,t); 

^/'(x, t) is a complex, classical field and g is the effective interaction strength [551156] 
This equation conserves the total energy 



E 



and the total number of particles 



A 



iv^r+-^?i^r 



(2) 



N= I \ip\'^d'^x 



(3) 



A 



where ^ = is the area of our 2D, periodic, computational domain of side L. From 
([1]) we obtain the continuity equation 



dp 
dt 



+ V ■ (pv) = 0, 



(4) 



where p = \if)\'^ is interpreted as the particle density and the velocity is 



(5) 
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We can use the Madelung transformation ip{yi,t) = y/pe*^^'''*^ where ^(x, t) is the phase 
of ■ip{'x,t), to write v(x, t) = 2V6'(x, t), whence we get |35] 



E 



A L 



d X — Ekin + Eint + Eq^ 



(6) 



where the kinetic, interaction, and quantum-pressure energies are defined, respectively, 
as 



4 J A 

V\p'/^\^d'x. 



(7a) 
(7b) 
(7c) 



A 



We separate the compressible (supercript c) and the incompressible (superscript i) parts 
of the kinetic energy by making use of the decomposition 



where V ■ (p^^^v)* = and V x (p^/^v)"^ = 0, whence we obtain the following: 

EL = ll \iVPv)fd\ 



El.n = ^jJi^/~P^)rd'x. 



The spectra for these energies are defined as follows: 

1 



E 



kin 



Ekin 



|(pi/2v)fd2A;= / El,^{k)dk; 



int 



|(pV2v)frf2A;^ j EUk)dk- 
\^/^^^\H^k = [ EUk)dk- 



and 



furthermore, we define an occupation-number spectrum n{k) via 



Eq = /|VpV2|2rf2^= Eg{k)dk; 



N 



\iP\'d'k = / n{k)dk; 



(8) 

(9a) 
(9b) 

(10) 

(11) 
(12) 

(13) 
(14) 



here we denote the Fourier transform of ^(x) by A; and, for notational convenience, we 
do not show explicitly the dependence of these spectra on time t. In any computational 
study, we must limit the number of Fourier modes that we use in our study of the GP 
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equation; we refer to such a GP equation as a Fourier-truncated GP equation (cf. |15|ll^ 
for studies of the Fourier- or Galer kin-truncated Euler equation). 

The Bogoluibov dispersion relation UB{k) is obtained by hnearizing ([T]) around a 
constant ip. For a total number of particles ([3]) = 1, it is 



u;B(k) = kc^l + ^, (15) 
where the sound velocity is c = and the coherence length is 

We investigate thermalization in the 2D GP equation, so it is useful to recall that 
a uniform, interacting, 2D Bose gas has a high-temperature disordered phase and a 
low-temperature, Berezenskii-Kosterlitz-Thouless (BKT) phase [571460] . which shows 
quasi-long-range order with an algebraic decay of the spatial correlation function [52] 

c(r) = ([e"''^W - (e"'^W)] [e'^(^+^^ - (e'''("+'-)))]); (17) 

for temperatures T below the transition temperature Tbkt (or energy Ebkt in the 
microcanonical ensemble), 

c(r) ~ r"'', (18) 

where r = |r| and the critical exponent t] < 0.25 for T < Tbkt] and rj = 0.25 at 
T = Tbkt [S3]- The BKT phase shows bound vortex-antivortex pairs; these unbind 
above Tbkt, so 

c(r) ~ e-'■/^ (19) 
in the disordered phase, with i the correlation length. 

2.2. Numerical Methods and Initial Conditions 

To perform a systematic, pseudospectral, direct numerical simulation (DNS) of the 
spatiotemporal evolution of the 2D, Fourier-truncated, GP equation, we have developed 
a parallel, MPI code in which we discretize ^(x, t) on a square simulation domain of side 
L = 32 with collocation points. We use periodic boundary conditions in both spatial 
directions, because we study homogeneous, isotropic turbulence in this 2D system, and 
a fourth-order, Runge-Kutta scheme, with time step At, for time marching. We evaluate 
the linear term in ([T]) in Fourier space and the nonlinear term in physical space; for the 
Fourier-transform operations we use the FFTW library [61]. Thus, the maximum wave 
number kmax = {Nc/'2)Ak, where Ak = 27r/L, and 

^kmax = (20) 
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We have checked that, for the quantities we calculate, dealiasing of our pseudospectral 
code does not change our results substantially; here we present the results from our 
pseudospectral simulations that do not use dealiasing. 

To initiate turbulence in the 2D GP equation we use three types of initial conditions 
ICl [H], IC2, and ICS [51], always normalized to correspond to a total number of 
particles ([3]) = 1. The first of these is best represented in Fourier space as follows: 

^(k,t = 0) = -^l=exp (^-fc-^^ expmK,ky)), (21) 

where k = ^/k^'+k^, Qik^^ky) are random numbers distributed uniformly on the 
interval [O, 27r] ; k^ = J\foAk and a = BAk, where the integer Ao controls the spatial 
scale at which energy is injected into the system, and the real number B specifies the 
Fourier-space width of ip at time t = 0. The initial condition IC2 is like ICl but, in 
addition, it has a finite initial condensate population A^^q =| ?/'(k = 0,t) p (Afc)^ at time 
t = 0. 

We obtain the initial condition ICS by solving the 2D, stochastic, Ginzburg-Landau 
equation (SGLE), which follows from the free-energy functional 

^ = j/^ (iV^P - ^I^P + \9\^\'^ , (22) 
where n is the chemical potentia|§|. The SGLE is 

where C is a zero-mean, Gaussian white noise with 

(C(x,t)C(x',t')) = ^5(x-x')5(t-0, (24) 

where D = 2T, in accordance with the fluctuation-dissipation theorem [62], T is the 
temperature, and 6 the Dirac delta function. Finally, the SGLE (|23l) becomes 

^ = V'^-^ij + g\^\^^ + C, (25) 

which we solve along with the following, ad-hoc equation 

f^-^(iV-A'„). (26) 

to control the number of particles A^; the parameter A'av controls the mean value of N; 
and governs the rate at which the SGLE equilibrates. We solve the SGLE by using 
a pseudospectral method, similar to the one described above for the 2D, GP equation, 
with periodic boundary conditions in space, an implicit-Euler scheme, with time step 
At, for time marching and the method of reference [63] (see page 25 of this reference). 

§ Recall that the SGLE can be thought of as an imaginary-time GP equation with external, additive 
noise (see, e.g. reference [38] ) 
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Table 1. Parameters for our DNS runs A1-A13, B1-B2, and C1-C6: N'^ is the 
number of collocation points, /cq is the energy-injection scale, a is the Fourier- 
space width of "0 at t = 0; 5 is the effective interaction strength; Nq is the 
initial condensate population; D and are respectively, the variance of the 
white-noise and the initial value of the truncation wave number, which we use 
in the initial conditions of type IC3; E is the total energy; we use a square 
simulation domain of area A = L"^] we choose L = 32. 



3. Results 

We first present the time evolution of the different energies, the probability distribution 
functions (PDFs) of the velocity components, and the population Nq in the zero-wave- 
number mode. We then give a detailed statistical characterization of the temporal 
evolution of the Fourier-truncated, 2D, GP equation in the four regimes mentioned in 
the Introduction (section [T]) . 

3.1. Evolution of energies, velocity PDFs, and the zero-wave-number population 

We show the early stages of the time evolution of the energies E\^^, -^fcm' -^mt, and Eq, 
from our DNS runs A1-A4, Bl, and C6 in figured] The runs A1-A4 use initial conditions 
of type ICl, in which E\^^ is a significant fraction of the total initial energy; the runs Bl 
and C6 start with initial configurations of type IC2 and ICS, respectively, in which E].^^ 
is negligibly small at t = 0. The transient nature of the early stages of the dynamical 
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Figure 1. Plots versus time t of the four components of the total energy E^-^, 
E^-^, Eint, and Eg, during the initial stages of evolution, from our DNS runs 
(a) Al, (b) A2 (c) A3, (d) A4, (e) Bl, and (f) C6 (see tabled]). 



evolution of the dissipationless, unforced, 2D, GP equation is evident from figure [H in 
which we observe a rapid conversion of El^^ into the other three components, with a 
significant fraction being transferred to E^-^^; moreover, the transient stage depends on 
the initial conditions, as we describe below. Figures [1] (a)- (c), show comparisons of the 
temporal evolution of the energies, from the runs A1-A3; we observe, in particular, that 
the conversion of El-^ into the other energy components is accelerated as g increases 
from 1000 to 5000 (cf. [12]); and there is a corresponding acceleration in the approach 
to thermalization. Moreover, the larger the value of El.^^ the larger is the time required 
for thermalization, as we can see by comparing figures [1] (a) and (d), for the runs Al 
and A4, respectively; the run A4 starts with a high value of El-^{t = 0) because of a 
large number of vortices and anti-vortices, so it takes a long time to thermalize; indeed, 
if the spatial resolution of our DNS is very high, the computational cost of achieving a 
statistically steady state is prohibitively high for intial conditions A1-A4. In contrast, 
the runs Bl and C6 have negligibly small values of El-^(t = 0) to begin with (figures [1] 
(e) and (f), respectively); and E\^^{t) remains close to zero throughout the dynamical 
evolution here. For run Bl, both E'^^^ and Eg start from values close to zero, grow at 
the cost of Einti and finally saturate to small, statistically steady values. For run C6, 
there are hardly any vortices in the initial configuration, so the energies start fluctuating 
about their statistically steady values very rapidly. 

In figure [2] we plot, at three instants of time, the PDFs of Vx and Vy, the Cartesian 
components of the velocity, for our DNS runs Al, Bl, and C6, which correspond, 
respectively, to initial conditions of types ICl, IC2, and ICS. For the run Al, these 
PDFs, in figures |2] (a)-(c), show a crossover from a distribution with power-law tails 
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to one that is Gaussian; the right and left tails of the PDFs in figure [2] (a) can be 
fit to the form ~ v~'^, with 7 ~ 3.2, and i = x 01 y (we show fits only for i = x). 
Such power-law tails in velocity-component PDFs have been seen in experiments [6l] 
and some numerical studies [391 HH [651 [66]. However, it has not been noted hitherto 
that, for turbulence in the Fourier-truncated, 2D, GP equation with low-energy initial 
conditions, such PDFs evolve, as t increases, from PDFs with power-law tails (figure [2] 
(a) for run Al), to ones with a Gaussian form near the mean, followed by broad tails 
(figure [2] (b) for run Al), and then to more-or-less Gaussian PDFs (figure [2] (c) for 
run Al), but with tails that can be fit to an exponential form. This evolution towards 
Gaussian PDFs is associated with the annihilation of vortices and anti-vortices. The 
Video SI in the Supplementary Material shows the temporal evolution of this PDF in 
the left panel and the spatiotemporal evolution of the pseudocolor plot of the vorticity 
in the right panel. The analogues of figures [2|(a)-(c) for runs Bl and CI, both of which 
have a negligibly small value of E^^^ at t = 0, are given, respectively, in figures [2]^d)-(f) 
and figures [21(g)- (i). 

We turn now to the time evolution of the population No{t), in the k = 
mode [36l[ini[67j, and its dependence on the initial conditions. In figure [3] (a) we plot 
A'o versus t for the runs A1-A4 (red, blue, green, and brown curves, respectively), which 
use initial configurations of type ICl; these figures show that N^it) increases with t, on 
average, and depends on E, g, ko, and a. For the runs Al and A2 (red and blue curves 
in figure [3] (a)), No(t) approaches a saturation value for the time scales probed by our 
simulations; figure [3] (a) also shows that, as we increase g (red, blue, and green lines in 
figure [3] (a)), the fluctuations in Nq are enhanced and its large-t value, which it seems 
to approach asymptotically, diminishes. By comparing the runs Al and A4 (red and 
brown lines in figures [3] (a)), we see that the latter has a higher value of E than the 
former, because both k^ and a are smaller for Al than for A4; thus, A'o(t) grows more 
slowly in A4 than in Al; and, after an equal amount of simulation time, its value in A4 
is nearly an order of magnitude lower than in Al; the former shows large fluctuations in 
No{t) and no sign of saturation. The run Bl (figure [3] (e)) uses an initial configuration 
of type IC2, with a large value of iVo(t = 0) = 0.95; in this case, after a period of initial 
transients, A^o(^) 0.98 over our simulation time. The run C6 (figure [3] (f)) uses an 
initial condition of type ICS; here No(t) fluctuates slightly but remains close to its initial 
value (cf. 

To study the dependence of A^o(^) on the number of collocation points A'^^, we evolve 
the initial configuration of Al for Nc = 512 (run A5), Nc = 256 (run A6), Nc = 128 (run 
A7), and Nc = 64 (run A8). Figure [3] (g) shows plots of No{t) versus t for these five runs; 
clearly, the initial evolution of No(t) depends significantly on A^^^; however, the large-t 
values of iVo(t), on the time scales of our runs, are comparable (^ 0.9) for the runs 
wth Nc = 128 (run A7), N^ = 256 (run A6), and N^ = 1024 (run Al). In contrast, the 
saturation value for the run with N^ = 64 (run A8) is ~ 0.8. For the run A5 {N^ = 512), 
No{t) shows large fluctuations and no sign of saturation over the time scale that we have 
covered; this suggests that Noit) also depends on the realisation of the random phases 
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Figure 2. Semilog (base 10) plots of the PDFs of the x (red circles) and y (green 
squares) components of the velocity from our DNS runs: (a)-(c) Al, (d)-(f) Bl, 
and (g)-(i) C6, corresponding to each of the three types of initial conditions ICl, 
IC2, and ICS, respectively. The complete time evolution of the PDFs (a)-(c) 
for the run Al is illustrated in the top- left panel of the video SI (supplementary 
material). The blue-dashed lines in (b)-(i) indicate fits to Gaussian PDFs; the 
dashed lines in (a) indicate power-law fits to the left (blue-dashed line) and 
right (orange-dashed line) tails of the PDFs (see text). 



Q{kx, ky) in f l2T]) . These plots of iVo(t) illustrate that complete thermalization proceeds 
very slowly for Nq\ in the completely thermalized state of the Fourier-truncated, 2D, 
GP system, A'^o must vanish in the thermodynamic limit by virtue of the Hohenberg- 
Mermin- Wagner theorem [SZllSB]; however, it is not easy to realize this limit in finite-size 
systems and with the limited run times that are dictated by computational resources. 
We discuss these issues again in section 13.41 and also refer the reader to [671|68] . 

3.2. Initial transients and the onset of thermalization 

The initial stages of the evolution of energy spectra for the Fourier-truncated, 2D, GP 
equations are qualitatively different for initial conditions of types ICl, 1C2, and 1C3. 
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t t 

Figure 3. Plots versus time t of the population Nq, in the zero- wave- number 
mode, from our DNS runs (a) A1-A4 (initial condition of type ICl), (b) Bl 
(initial condition of type IC2), (c) C6 (initial condition of type IC3), and (d) 
Al and A5-A8, for five values for the number of collocation points N^, namely, 
10242, 512^, 256^, 128^, and 642. 

The first type begins with a sizeable incompressible kinetic energy spectrum E\^^{k)] 
and the initial transients are associated with the annihilation and creation of vortex- 
antivortex pairs, the associated depletion of £'^^^(/c), and the growth of the other energy 
components [H]. In contrast, runs with initial conditions of types IC2 and ICS start 
with a very small incompressible-energy component, therefore, even the early stages of 
their dynamical evolution are akin to the late stages of the dynamical evolution with 
initial conditions of type ICl. In figures H] (a)-(d) we show the time evolution of the 
spectra E\^^{k)^ for the runs Al, A2, A3, and A4, to ascertain the presence of scaling 
behaviour, if any. We find that, in the low-Zc region, E]^^^{k) lacks a well-defined scaling 
region (unlike in |l2]); indeed, this region depends on the initial configuration, changes 
continuously with time, and, in particular, a k~^/'^ scaling region is tenable (a) over a 
range of wave numbers that is very tiny and (b) over a fieetingly short interval of time 
(around t = 50 for the run Al). At large wave numbers, El-^{k) ~ k~^, during the 
initial stages of evolution, because of the presence of the vortices this power-law 

form holds over the same time scales for which the PDF P{vx/<y^^) ~ v~"' (figures [2] 
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Figure 4. Log- log (base 10) plots of the spectra El.-^{k) from our DNS runs 
(a) Al, (b) A2, (c) A3, and (d) A4 at different times t (indicated by curves 
of different colours); a A;~^ power law is shown by orange-dashed lines. The 
complete time evolution of the spectra in (a), (b), (c), and (d) is illustrated in 
the video S2 (supplementary material). 



(a)-(b)). 

The initial transients described above are followed by a regime in which the energy 
and occupation-number spectra begin to show power-law-scaling behaviours, but the 
power-law exponent and the extent of the scaling region change with time and depend 
on the initial conditions; we regard this as the onset of thermalization, which is shown 
in figures |5] and [6l where we illustrate the time evolution of E^^^. Figure [5] (a) shows 
E^-^{k) for the run Al; we begin to see a power-law region here with Ej^^^{k) ~ k, on 
the \ow-k side of the peak after which the spectrum falls steeply. A similar E^^^ik) ~ k 
behaviour starts to emerge in the region k < k^ax for the run Bl (figure E] (g)). In this 
onset-of-thermalization regime, we also see the development of the following power laws: 
Eint{k) + Eq{k) ~ k (figure E]) and n{k) ~ l/k (figure E]). 

3.3. Partial thermalization and self -truncation 

3.3.1. Partial thermalization In the third stage of the dynamical evolution of the 2D, 
Fourier-truncated, GP equation, which we refer to as the partial-thermalization stage. 
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Figure 5. Log-log (base 10) plots of the spectra E^^^{k) from our DNS runs 
(a)-(c) Al, (d)-(f) A4, and (g)-(i) Bl at different times t (indicated by curves of 
different colours); a k power law is shown by orange-dashed lines. 

well-defined, power-law-scaling behaviours appear in energy and occupation-number 
spectra, with exponents that are independent of the initial conditions as illustrated by 
the compressible- kinetic-energy spectra in figures O (b), (c) (e), (f), and|n](b) for initial 
conditions of type ICl, and figures [5] (h) and[6] (e), and (f), for initial conditions of type 
IC2. It is important to distinguish between (I) spectra that fall steeply at large values 
of k, e.g., the spectra in figures [S](b), (c) (e), (f), and|n](e) and (f), and (ll) spectra that 
increase all the way to kmax, e.g., the spectra in figures |5] (h) and[6](b) and (c). In case 
(I), we have spectral convergence to the 2D GP partial differential equation (PDE); in 
case (II), the effects of Fourier truncation are so pronounced that our truncated 2D, GP 
system does not provide a good representation of the 2D, GP PDE. As we show below, 
case (I) can be further subdivided into (A) a subclass in which the maximum, at k — k^ 
in Ekin{k) = (Ej^ini^) + ^lini^))^ referred to as the self-truncation wave number [5Tj . 
moves out to kmax as a power of t and (B) a subclass in which k^ m.oves out to k^ax at 
a rate that is slower than a power of t. 

Figures (g)-(i), from the run Bl, show how E^^^ik) evolves as the spectral 




Figure 6. Log-log (base 10) plots of the spectra E^^^i^) our DNS runs 

(a)-(c) A7 and (d)-(f) B2 at different times t (indicated by curves of different 
colours); a k power law is shown by orange-dashed lines. 

convergence to the GP PDE is lost in case (II); note that the scahng region with 
^kin ~ ^ s^^s ^is'^ wave numbers close to kmax and then extends to the low-wave- 
number regime. For case (lA) analogous plots of El^^{k) are given in, e.g., figures [6] 
(a)-(c). We give plots for case (IB) in the next subsection, where we study in detail the 
time dependence of kc- Illustrative plots of the spectra {Ei{k) + Eq{k)) and n{k) in this 
regime of partial thermalization are given in figures [7] and IH respectively. 

3.3.2. Self -truncation We now present a detailed characterization of the partial- 
thermalization regime, when energy spectra display self-truncation at wave-numbers 
beyond fcc(i), which can be defined as follows: 

as the system approaches complete thermalization, kc{t) — )■ k^ax- In particular, we 
explore how the scaling ranges in energy spectra grow with t for different values of g, 
with the initial configuration and number of collocation points held fixed. For an 
initial condition of type ICl, with ko = 5Ak, a = 2Ak, and = 256, we obtain the 
time evolution of energy spectra for g = 1000 (run A6), g = 2000 (run A9), and g = 5000 
(run AlO) in figures [9] (a), (b), and (c), respectively, and their video analogues (Videos 
S3 (panel V2) in the Supplementary Material). The larger the value of g, the more 
rapid is the thermalization, and the consequent loss of spectral convergence, as we can 
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Figure 7. Log-log (base 10) plots of the spectra Eint{k) + Eq{k) from our DNS 
runs (a)-(c) Al, (d)-(f) A4, and (g)-(i) Bl at different times t (indicated by 
curves of different colours); a k power law is shown by orange-dashed lines. 

see by comparing the sky-blue (run AlO), green (run A9), and purple (run A6) spectra 
in figures |9] (a)- (c); run A6 loses spectral convergence around t = 2500. We obtain the 
same qualitative g dependence, with /cq = 15A/c, a = 2Ak, and = 256, for g = 1000, 
2000, and 5000, i.e., runs All, A12, and A13, respectively, for which energy spectra are 
portrayed in figures [9] (d)-(f) and Video S3 (panel V3) in the Supplementary Material. 

In figures [3 (g)-(i) we explore the Nc dependence of the self-truncation of energy 
spectra, for initial conditions, with ko = 5Ak, a = 2Ak, and g = 1000, and five different 
values of N^, namely, A^^ = 1024 (run Al), 512(run A5), 256 (run A6), 128 (run A7), and 
64 (run A8). We find, not surprisingly, that the lower the value of the more rapidly 
does the system lose spectral convergence. 

Initial conditions of type IC2 lead to energy spectra whose time evolution, and their 
dependence on g and Nc, is similar to those that are obtained from intial conditions of 
type ICl. 

With initial conditions of types ICl and IC2, we cannot control the initial value 
kc{t = 0) = /c*" easily. However, initial conditions of type ICS, which we obtain from 
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Figure 8. Log-log (base 10) plots of the spectra n{k) from our DNS runs (a)- 
(c) Al, (d)-(f) A4, and (g)-(i) Bl at different times t (indicated by curves of 
different colours); a. k^^ power law is shown by orange-dashed lines. 



the SGLE, allow us to control /c*" and start, therefore, with initial spectra that display 
partial thermalization for k < fc*" [51] and a sharp fall thereafter. In figure [TO] we show 
the time evolution of E'j:-^{k) for such initial conditions from runs C1-C6. For different 
representative values of /c*", g, and D, we now study the time evolution of kc{t), which 
characterizes the growth of the partially thermalized scaling region. Here too, as with 
initial conditions of types I CI and IC2, if all other parameters like fc*" = 6.0 and D are 
held fixed, the speed of thermalization increases with g (cf. figure [10] (a) for the run CI, 
with g = 5000, and figure [TO] (b) for the run C2, with g = 1000). For these runs C1-C6, 
the growth of the energy spectra, in the region k > /c*", starts with the smoothening 
of the sharp cut-off at /c*"; the higher the value of /c*", the slower is this growth (cf. 
figures [To] (b), (d), (e), and (f) for runs C2, C4, C5, and C6, respectively). By contrast, 
an increase in D (or T) in the SGLE, accelerates this growth (cf. figures [TO] (b) and (c) 
for runs C2 and C3, respectively). 

The growth of kc{t) with t, illustrated in figure [TT] (a), can be fit to the form 
kc{t) ~ t"; however, as we show below, a depends on the initial condition. We obtain 
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Figure 9. Log-log (base 10) plots of the spectra Ekinik) from our DNS runs 
(a)-(c) A6, A9, and AlO (fco = 5AA; and a = 2Ak), (d)-(f) All, A12 and A13 
(ko = 15Ak and a = 2A), and (g)-(i) Al, A5-A8 (N^ = 1024^, 512^, 256^, 128^, 
and 64^). The complete time evolutions of the spectra in (a)-(c), (d)-(f), and, 
(g)-(i) are illustrated in the panels V2, V3, and V4 of video S3. 



the exponent a either from slopes of log-log plots of (i) kc{t) versus t or (ii) dkc/dt versus 
kc/kmax] we denote the values from procedures (i) and (ii) as ai and ^2, respectively. 
Note that in (ii) we have a parametric plot [38l|51], shown in figure [11] (b); this yields a 
straight-line scaling regime with slope x 02 = 1/(1 — x)- The values of ai and 0:2, 
listed in tabled show that ai ~ 0:2; the discrepancy between these two values for a is 
a convenient measure of the errors of our estimates. For runs C4, C5, and C6, we cannot 
obtain 02 reliably; the small values of ai for these runs indicate very slow growth of 
kc(t); indeed, in runs C5 and C6, a case can be made for a logarithmic growth of kc(t) 
with t. 

3.4. Complete thermalization 

The partially thermalized stage of the dynamical evolution of the 2D, Fourier-truncated, 
GP equation may either gradually become completely thermalized, in which state a 
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Figure 10. Log- log (base 10) plots of the spectra E]^in{k) from our DNS runs 
(initial conditions of type ICS (a) CI, (b) C2, (c) C3, (d) C4, (e) C5, and (f) C6. 




Figure 11. Plots of (a) the self-truncation wave- number kc{t) versus time t and 
(b) dkc/dt versus kc/kmax from our DNS runs A1-A4, B2, and C1-C6. 



power-law scaling region is present in the entire energy and the occupation number 
spectra, or remain self-truncated with logarithmic growth. In figures (g)-(i) and E] 
(a)-(c), we show the compressible kinetic energy spectra E^^^ for the runs Bl and A7, 
where E^^^ shows power-law scaling over the entire wave number range, from k = 27t/L 
up to kmax, towards the end of the respective simulations; a naive fit is consistent with 
E^-^{k) ~ k (but see below). 

3.4-1 ■ Correlation functions and the BKT transition A uniform, 2D, interacting Bose 
gas exhibits a BKT phase at low energies (temperatures in the canonical ensemble). 
Thus, the completely thermalized state of the 2D, Fourier-truncated, GP equation 
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Table 2. Summary of the self-truncation results from our DNS runs A1-A4, B2, 
and C1-C6: E is the total energy; kmax = 2'nNc/2L; ^ = L/y/g is the healing 
length; and kl are the initial and final values of kc (averaged over a few time 
steps); ai is the slope obtained from the log- log (base 10) plot of kc versus t 
and a2 = 1/(1 — x); where x is the slope obtained from the log-log (base 10) 
plot of dkc/dt versus kc/k^ax- 

should yield a BKT phase [52ll5l] . with the correlation function c(r) ~ r"'', at energies 
E < Ebkt', and c(r) should decay exponentially with r if E > Ebkt- We show this 
explicitly now by using initial conditions of type ICl with Nc = 64 and Nc = 128 and 
g = 1000; we obtain different energies by changing kg and a (runs D1-D13 and E1-E12 
in table |3]). 

In figure[T2l we present plots of the correlation functions c(r). To illustrate the BKT 
transition clearly, we present log-log plots of c(r) versus r, for E < -Ebkt; in figures [12] 
(a) and (d), where the straight lines indicate power-law regimes; and, for E > Ebkt, 
we use semi-log plots, as in figures [12] (b) and (e), where the straight lines signify an 
exponential decay of c(r) with r. Given the resolution of our DNS runs, we find that, 
in a small energy range in the vicinity of -Ebkt, we cannot fit power-law or exponential 
forms satisfactorily; this leads to an uncertainty in our estimate for -Ebkt- Aside from 
this uncertainty, the behavior of c(r), in the regime of complete thermalization, is in 
accord with our expectations for the BKT phase; in particular, the exponent r] (see 
equation ( IT8l) ) depends on E for E < Ebkt as shown in figures [I2] (c) and (f). Our 
values for 77, for the runs with E < Ebkt and with = 64 and = 128, are listed 
in table [3] Note that Ebkt ^ Q{Nc = 64) and Ebkt ^ 19(A^c = 128), i.e., Ebkt 
depends on Nc, the number of collocation points; we show analytically below how a low- 
temperature analysis can be used to understand this dependence of Ebkt on Nc- In the 
completely thermalized state of the Fourier-truncated, 2D, GP system, A^o must vanish 
in the thermodynamic limit by virtue of the Hohenberg-Mermin- Wagner theorem [57i[58] 
and n{k) ~ k~^~^^; it is not easy to realize this limit in finite-size systems and with the 
limited run times that are dictated by computational resources (see the plots of Nq in 
figure [3]); however, finite-size scaling can be used to extract the exponent rj from the 
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Table 3. List of parameters for our complete-thermalization DNS runs D1-D13 
{N^ = 128^) and E1-E12 (A^^ = 542-). jy2 jg ^^le number of collocation points; 
ko is the energy-injection scale; a is Fourier-space width oi ijj at t = 0; E is the 
total energy; and ij is the exponent of the correlation function c(r) ~ for 
E < Ebkt- S = 1000 for all the DNS runs and they have been performed on 
a square simulation domain of area A = L^, with L = 32. 



k = part of n{k) as shown in reference [68]; similarly, E^^^{k) should also show a 
power-law form with an exponent that depends on 77, but this is difficult to realize in 
numerical calculations with limited spatial resolutions and run lengths. 

3.4-2. Analytical estimation of the energy of the BKT transition The energy of a pure 
condensate of a uniform, weakly interacting, 2D Bose gas, which is described by the GP 
equation ([1]), is Eq = g/ {2A). We define the energy of our system to be = Eq{1 + 6S); 
this energy E is fixed by the initial condition; and 6S measures the relative amount by 



which E exceeds Eq. As we show in the Appendix B the Nc dependence of the energy 



-E'bkt, at which the BKT transition occurs, can be obtained approximately as follows. 
We begin with 

8 

S^BKT = S^BKT 7 (^^ 

log(7r2iV,2 



29 



where SSbkt, the estimate for the BKT transition energy that follows from an energy- 
entropy argument (see fEUl) in the Appendix and is 

2 AT 2 f2p 

whence we obtain 

SSbkt = ^. (30) 
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Figure 12. Plots of c(r) versus r for different energies in the complete- 
thermalization regime, for = 128^ ((a) and (b)) and = 64^ ((d) and 
(e)). (a) and (d) Log-log (base 10) plots of c(r) versus r for different energies 
E < Ekt'i the slopes of the linear parts of these plots yield the exponent 
rj (table [3|); (b) and (e) semilog (base(lO) plots of c(r) versus r for different 
energies E > Ekt\ (c) {N"^ = 128^) and (f) [N'^ = 64^) show plots of and 
Nq versus E (on the time scales of our runs Nq is nonzero; see the text for a 
detailed discussion). 
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39.44 
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Table 4. The values of Eq, 5£bkt (see 5£ (see ([30])), £^^kt (see ([3T]) ). 

and from our DNS runs D1-D13 {N^ = 64) and E1-E12 {Nc = 64). Eq 

is the ground state energy of a pure condensate of a uniform, interacting, 2D 
Bose gas and -E'bkt BKT-transition energy determined using our DNS runs. 



We can now write 

by using this expression we can determine the ratio -E'bkt(^c)/-^bkt(^c) for runs with 
two different values, and iV^, for the number of collocation points; we can also 
obtain this ratio from our DNS, by determining the value of E at which the exponent 
7] becomes 1/4. In Table H] we compare E-QwriNc) for = 64 and Nc = 128; our 
analytical approximation (13T|) yields E^^r^./ E^rj, ^ 3.15; this is in excellent agreement 
with the value ~ 3.14 that we obtain for this ratio from our DNS results. 
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4. Conclusions 

We have carried out an extensive study of the statistical properties of the dissipationless, 
unforced, 2D, Fourier-truncated, GP equation. Our study has been designed specifically 
to study and identify the universal features, if any, of the turbulent evolution of the 
solutions of this equation, by undertaking a systematic DNS. In our study, we have 
used statistical measures such as velocity-component PDFs and energy and occupation- 
number spectra, for a large number of initial conditions. To the best of our knowledge, 
such a comprehensive study of the Fourier-truncated, 2D, GP equation has not been 
attempted hitherto. 

Our comprehensive study of the Fourier-truncated, 2D, GP equation, which makes 
use of the three types of initial conditions (section 12. 2p and a wide range of parameters 
(tables[T]and|3]), allows us to systematize the dynamical evolution of this system into four 
different regimes, with qualitatively different statistical properties. This demarkation 
of the evolution into different regimes has not been systematized in earlier studies, 
which have concentrated only on one or two of these regimes. For example, the study 
of reference [39] has investigated states with a significant number of vortex-anitvortex 
pairs and obtained for them PDFs of velocity components that have power-law tails of 
the type shown in figure [2j References [171[5lll68] have investigated the BKT nature 
of the thermalized state. Wave-turbulence studies [32],|1T],|69] have focussed on power- 
law regions in energy and occupation-number spectra of the type we find in our third 
regime. The DNS studies in |¥TI - H5l[70] have considered the time evolution of spectra 
and PDFs for the Fourier-truncated, 2D, GP equation; in some cases, these studies 
introduce dissipation or hyperviscosity and forcing; they have also reported different 
power laws in spectra |l2 |H3lH5] . Our work suggests that, at least in the dissipationless, 
unforced, Fourier-truncated, 2D, GP equation, the only robust power laws in spectra 
are the the ones we have reported above; all other apparent power laws occur either (a) 
for very special initial conditions [H] or (b) last for fleetingly small intervals of time and 
extend over very small ranges of k. 

To recapitulate, we find that, in the first dynamical-evolution regime of the Fourier- 
truncated, 2D, GP equation, there are initial-condition-dependent transients. In the 
second regime the energy and the occupation-number spectra start to develop power- 
law scaling regions, but the power-law exponent and the extent of the scaling region 
change with time and are influenced by the initial conditions. In the third regime, 
of partial thermalization, we find E1^^{k) and Eint{k) + Eq{k) ~ k, and n{k) ~ l//c, 
for k < kc{t) and, for k > kc, we find an initial-condition-dependent self-truncation 
regime, in which the spectra drop rapidly; the self-truncation wave number kc{t) grows 
either as t° or logarthimically for different intial conditions (table [2]). In the fourth, 
complete-thermalization regime, power-law forms of correlation functions and spectra, 
for E < -Ebkt; are consistent with their nontrivial BKT forms; however, considerable 
care must be exercised, as explained in section [3.4. II and [4711541168] . to distinguish these 
nontrivial power laws from their wave-turbulence analogs [32| 1 ^ 169]. 
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Appendix A. 

The GP equation, which describes the dynamical evolution of the wave function '?/'(x, t) 
of a weakly interacting, 2D Bose gas at low temperatures, is 

= -|^V^^(x,t) +^2D|V^P^(x,t), (A.l) 

where g2D is the effective interaction strength. As we have mentioned earlier (see ([2]) 
and ([3])), the GP equation conserves the energy, given by the Hamiltonian 

and the total number of particles n = {ipl'^d'^x. To obtain ([T]), we first divide (lA.ip 
by h and define g = 5'2d/^; we then set h/2m = 1, with m = 1, so that l^/'p is the same 
as p; this is tantamount to using units with h = 2. 



Appendix B. 

The Berezinskii-Kosterlitz-Thouless (BKT) transition is best studied by using the 
renormalization group |52]; here, we restrict ourselves to the heuristic, energy-entropy 
argument to obtain a rough estimate of the BKT transition temperature Tbkt- In the 
XY model, this transition is studied by using the Hamiltonian 

^xy = -J5^cos(^,-%), (B.l) 

<i,j> 

where < i,j > denotes nearest-neighbour pairs of sites, on a 2D square lattice, J is the 
nearest-neighbour exchange coupling, and {6i — 9j) is the angle between the nearest- 
neighbour, XY spins on sites i and j. In the continuum limit, the above Hamiltonian 
becomes, to lowest order in spatial gradients, 

HxY = ^l dM'^Oix))'. (B.2) 

By comparing ( IB. 21) with the kinetic-energy term in ( 1A.1I) . we find that 

m (27r) 

where P denotes the Onsager-Feynman quantum of velocity circulation P = 47r/i/2m 
h/m. A rough estimate for the BKT transition temperature Tbkt is given below: 

nJ Tx I (V^) |2 pp2 



■BKT 



2kB 2mkB Snks' 



(B.4) 
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here Tbkt denotes the estimate for Tbkt that follows from an energy-entropy 
argument [52]. For T < Tbkt? the phase correlation function c(r) (see f|T7|) ) and the 
angle- integrated spectrum c{k), which follows from a Fourier tranform of c(r), scale as 

T 

c(r) ~ (a/r) ■*'^BKT (B.5) 

and 

c{k)^k^^^^^, (B.6) 

respectively. Above Tbkt the correlation length 

Jk-^E{k)dk 
~ jE{k)dk ^ 

is finite; and, as T — )■ Tbkt, it displays the essential singularity 

I ~ exp(6(TBKT/(T - Tbkt))'/'). (B.8) 

Appendix B.l. 

We now develop an analytical framework, which is valid at low-temperatures T ^ Tbkt, 
that can be used to test some of the results of our DNS runs in the region of 
complete thermalization. We first calculate equilibrium thermodynamic functions for 
a weakly-interacting, 2D Bose gas, in the grand-canonical ensemble; we then obtain 
their analogues in the microcanonical ensemble. In the grand-canonical ensemble the 
probability of a given state is 

P=le-/3(^^-M^), (B.9) 

where S is the grand partition function, /3 the inverse temperature, /i the chemical 
potential, and the number of bosons. The grand-canonical potential is 

= -/3-Mog(S); (B.IO) 

and the mean energy E, entropy S, and are 

dVL 

N=--, (B.lla) 

S = p^dVL/dp, (B.llb) 



We adapt to 2D the 3D study of Ref. [51], expand i}) in terms of Fourier modes A^, and 
obtain VL as the sum of the saddle-point part VLgp and fig, the deviations from the saddle 
point that are quadratic in A^. We write Vt = Q^p + ^q, where Qgp = —Afi^/2g and 
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We can also calculate the condensate depletion 6N, where the particle number 
Nq + 6N and Nq is the number of particles in the k = mode, as follows: 



The integrals in the flB.12p and flB.lSp can be performed analytically, but, in contrast 
to the 3D case where the primitives are zero at p = 0, the 2D primitive for Vt^^ is finite 
at p = and for 5N it is infra-red (I.R.) divergent. By subtracting the I.R. finite and 
divergent terms from VLq and 5N , respectively, we get the following expressions, in 2D, 
in the thermodynamic limit oc: 



n 



2g A7r(3ff 27r(3K^ 
and 

mA (log(l + ^) + log^^ 
By using the thermodynamic relations (IB. lip , we get 



(B.14) 



r'max-^ 



= ^ jz^, (B.15) 



.2 ^ 
max 



g 27T(3h'^ 



N=^ °- ^'"M- (B.16) 



and 



^=^ + 4^^ 2^^ • ^^-^^^ 



Appendix B.2. 



We next determine the chemical potential /i, which fixes the total density p = niN/A 
at a given value, by solving the equation 



2 

mil log(l + T^iax) 

P- — + =0; (B.18) 



g 2'Kpn 
at /9 = oo, i.e., zero temperature (subscript 0) we obtain 



to order /3 we get 
where 



/xo = — ; (B.19) 

m 



fi = fio + Sfi, (B.20) 



mg {Agp^ + ppl^^^) log(l + 



'max 
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We insert /i from (lB.20p into (lB.15p . define the change in density Sp = mSN/A, use the 
energy E from flB.171) . and then expand to order to obtain 



m2 (^log(l + ^) + log(^ 



= ^ — 

and ^ ^ 

By using (lB.4p and p = m \ (ip) we obtain 

/3bkt = = (B.24) 



which we can use along with (IB.22P to relate the condensate relative depletion Sp/p to 



(3/(3bkt, where (3 = 1/(/cbT) and k-Q is the Boltzmann constant, as given below: 



We use this low-temperature result flB.25p to estimate the inverse-temperature scale 
/3bkt, at which the depletion of the k = condensate mode becomes significant for 
a finite-size system with A'^^ collocation points (which fixes the maximum momentum 
Praax)] in particular, we can solve ( ]B.25p . for Sp/p = 1, to obtain 



5bkT 1 /Pma.(l + ^)^\ 

By making the replacements that correspond do defining h, m, and g in terms of c and 
^, as in [51], Pmax ^^max, ^ ~> V^cm^, and g — c^m? / p, we can rewrite (1B.26P as 



/3bKT 1, A 2 /I /I , ^max^'C^A /o o'tN 

3 = o log ^max^ (1 + ^ ) • (B.27) 

Appendix B.3. 

Our DNS runs, which use initial conditions of type ICl and IC2, give the dynamical 
evolutions of the Fourier-truncated, 2D GP equation, which is a Hamiltonian system. 
The energy E, particle number iV, and area A are conserved in this evolution, so 
our calculation can be viewed as a simulation of this Hamiltonian system in the 
microcanonical ensemble, which yields, eventually, the fully thermalized state that we 
have described above. Therefore, we now transform the results, which we have obtained 
in the previous subsection, into their counterparts in the microcanonical ensemble. In 
the low-temperature limit, (IB.23P yields 

Iirh^ [Zm^E — gp A) 
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The energy of a pure condensate is 



/3-5.00 2 



Eo = lim E = (B.29) 



and the energy and the inverse temperature /3 ( 1B.28I) can be related as follows: 



E = Eo{l + S£), (B.30) 
where 6S is the relative increase of energy above Eq, and 

2 2 

If we now substitute /3 = (3bkt by using (lB.26p . we obtain, in terms of c, ^ and p (see 
text just below fiRM ) 

i^o = (B.32) 
5^BKT = (B.33) 



and 



5^BKT = 7 '"'^"^ (B.34) 



log (^lax^(l + %^)) ■ 



All the energies mentioned in the main paper are dimensionless; thus, to convert the 
energies given in this Appendix to dimensionless forms, we divide them by h. Hence, the 
energy of a pure condensate is obtained, in the dimensionless form, by dividing flB.29p 
by h, which gives 
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